function [IA_Network, consistent, iterations] = IA_Network_Path_Consistency_Hogge(IA_Network, i, j, IA_Relation)
%IA_Network_Path_Consistency To Add R(i, j)
%   Detailed explanation goes here

% Todo (RF): both 'N' and 'R' matrices are at least 2x bigger than required 
% as the relationship matrix is symmetric around the (top left to bottom 
% right) diagonal. The diagonal also contains no information as each 
% interval can only be equal to itself and no other relation is possible.

consistent = 1;
iterations = 0;

if IA_Relation_Is_Dummy(IA_Relation)
  return;
end

if IA_Relation_Is_Null(IA_Relation)
  consistent = 0;
  return;
end

N = IA_Network.relations;
N(i, j) = IA_Relation_Intersect(N(i, j), IA_Relation);
N(j, i) = IA_Relation_Invert(N(i, j));

if IA_Relation_Is_Null(N(i, j))
  consistent = 0;
  return;
end

% Add <i,j> to queue ToDo;
ToDo = struct('i', i, 'j', j);

% While ToDo is not empty do
% begin
while (size(ToDo, 2) > 0)
  % Get next <i,j> from queue ToDo;
  curr = ToDo(1);
  i = curr.i;
  j = curr.j;
  ToDo = ToDo(2:size(ToDo, 2));
  HasGreaterIJ = IA_Relation_Has_Greater(N(i, j));
  HasLessIJ = IA_Relation_Has_Less(N(i, j));
  HasContainsIJ = IA_Relation_Has_Contains(N(i, j));
  HasDuringIJ = IA_Relation_Has_During(N(i, j));
  % N(i, j) <- R(i, j)
  % N(j, i) <- R(j, i)
  % For each node k such that Comparable(k, j) do
  % begin
  for k = 1:size(IA_Network.nodes, 2)
    if IA_Interval_Comparable(IA_Network.nodes(k), IA_Network.nodes(j))...
    && (k ~= i) && (k ~= j)
      if (IA_Relation_Has_Less(N(k, i)) && HasGreaterIJ)...
      || (IA_Relation_Has_Greater(N(k, i)) && HasLessIJ)...
      || (IA_Relation_Has_During(N(k, i)) && HasContainsIJ)
        continue;
      end
      % R(k,j) <- R(k,j) I N(k,j) I Constraints(N(k,i), R(i,j))
      temp = IA_Relation_Intersect(N(k, j),... 
          IA_Relation_Constraints_Hogge(N(k, i), N(i, j), N(k, j)));
      % Check consistency  
      if IA_Relation_Is_Null(temp)
        consistent = 0;
        return;
      end
      % if R(k,j) C N(k,j)
      if IA_Relation_Proper_Subset(temp, N(k, j))
        % then add <k, j> to ToDo;
        ToDo(size(ToDo, 2) + 1) = struct('i', k, 'j', j);
        N(k, j) = temp;
        N(j, k) = IA_Relation_Invert(N(k, j));
      end
    end
  end
  % end
  % For each node k such that Comparable(i, k) do
  % begin
  for k = 1:size(IA_Network.nodes, 2)
    if IA_Interval_Comparable(IA_Network.nodes(i), IA_Network.nodes(k))...
    && (k ~= i) && (k ~= j)
      if (HasLessIJ && IA_Relation_Has_Greater(N(j, k)))...
      || (HasGreaterIJ && IA_Relation_Has_Less(N(j, k)))...
      || (HasDuringIJ && IA_Relation_Has_Contains(N(j, k)))
        continue;
      end
      % R(i,k) <- R(i,k) I N(i,k) I Constraints(R(i,j), N(j,k))
      temp = IA_Relation_Intersect(N(i, k),...
          IA_Relation_Constraints_Hogge(N(i, j), N(j, k), N(i, k)));
      % Check consistency
      if IA_Relation_Is_Null(temp)
        consistent = 0;
        return;
      end
      % if R(i,k) C N(i,k)
      if IA_Relation_Proper_Subset(temp, N(i, k))
        % then add <i, k> to ToDo;
        ToDo(size(ToDo, 2) + 1) = struct('i', i, 'j', k);
        N(i, k) = temp;
        N(k, i) = IA_Relation_Invert(N(i, k));
      end
    end
  end
  % end
  iterations = iterations + 1;
end
% end;

IA_Network.relations = N;
end